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ABSTRACT 

We report on a series of numerical simulations of gas clouds with self-gravity forming 
sink particles, adopting an isothermal equation of state to isolate the effects of gravity 
from thermal physics on the resulting sink mass distributions. Simulations starting 
with supersonic velocity fluctuations develop sink mass functions with a high-mass 
power-law tail dN/d\ogM oc , T = — 1 ± 0.1, independent of the initial Mach 
number of the velocity field. Similar results but with weaker statistical significance 
hold for a simulation starting with initial density fluctuations. This mass function 
power-law dependence agrees with the asymptotic limit found by Zinnecker assuming 
Bondi-Hoyle-Littleton (BHL) accretion, even though the mass accretion rates of indi¬ 
vidual sinks show significant departures from the predicted M oc behavior. While 
BHL accretion is not strictly applicable due to the complexity of the environment, we 
argue that the final mass functions are the result of a relative dependence result¬ 
ing from gravitationally-focused accretion. Our simulations may show the power-law 
mass function particularly clearly compared with others because our adoption of an 
isothermal equation of state limits the effects of thermal physics in producing a broad 
initial fragmentation spectrum; L —>■ — 1 is an asymptotic limit found only when sink 
masses grow well beyond their initial values. While we have purposely eliminated many 
additional physical processes (radiative transfer, feedback) which can affect the stel¬ 
lar mass function, our results emphasize the importance of gravitational focusing for 
massive star formation. 

Key words: stars: formation — stars: luminosity function, mass function — ISM: 
clouds 


1 INTRODUCTION 

The origin of the stellar initial mass function (IMF) has 
been the subject of many theoretical studies over the years, 
with little resulting consensus . As discussed in the recent 
review by lOffner et al.l LolJl . currently popular theories 
may be classified into two groups: one in which masses 
of dense structures or cores created as a result of (super- 
sonic) turbulence map into t he resulting stellar masse s (e.g. , 
Padoan fc Nordlundl I 2 OO 2 I : iHennebelle fc Chabrieil l20ok 
2OO9I I. while the other invokes accretion from an extended 
region with dynam ical interactions betw ee n self-gravitating 
gas and stars le.g.. iBonnell et 'aDkOOl J U: iBate et ahlkoPa : 
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IClark et al.l l2007f ). Although there seems to be general 
agreement that thermal physics plays an essential role in 
the turnover of the IMF in the solar-mass (and below) 


reeime (iJaDDsen et al. 

I2OO5I: IBonnell et al.ll2006l:lBatell2012l: 

iKrumholz et al.l 2010l: 

Mvers et al.ll2014), the application of 


Jeans mass-type arguments is far less clear for massive stars, 
where continuing accretio n from regions beyo nd the local 
Jeans lengths is plausible (iBonnell et al.ll2001all . 


One plausible explanation for the upper mass IMF is 
gravitationally-foc used mass addition, sometimes referred to 
as “Bondi-Hoyle” (IZinneckeil Il982l : iBonnell et al.l 2001al ^ 
or “Bondi-Hoyle-Littleton” accretion (BHL; lEdgarl 2004l l . 
with an accretion rate given by 
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AtvG^M^Po 

(c2+„2)3/2 


= aM^ , 


( 1 ) 


where M is the mass accretion rate, M is the gravitating cen¬ 
tral mass, and po, Cs, and Vaa are the ambient gas density, 
the (isothermal) sound spee d, and the relativ e velocity of 
the ambient gas, respectively. IZinneck^ Jigsd i showed ana¬ 
lytically that, if a is constant, a population with a small but 
finite difference in initial masses will evolve to an asymptotic 
power law slope of F = dlog N/d log M = — 1. This is a sug¬ 
gestive result, given the apparent universality (within obser¬ 
vational uncertai nties) of the Salpete r slope F^ ~ —1.35 for 
high-mass stars llBastian et al.l 120101 ). and recognizing that 
the F —>■ — 1 is an asympotic limit reac hed only when t he 
accreted masses far exceed initial values llZinneckeilll982l l. 

However, the applicability of BHL accretion in mas- 
sive star formati o n has been challenged recently by 
Mjgchberger et ^ ll2014l . M14); (see also IPadoan et al.l 
2014 1. M14 analyzed the results of formation of sink par¬ 
ticles in a hydrodynamic simulation and argued that the 
accretion rates are generally M oc a result pre- 

dicted for the case of accretion in a gas-dominated potential 
llBonnell et ahlBoOlallO) rather than the the rate of equa¬ 
tion 0, w hich applies to accret ion in a “stellar-dominated 
potential” llBonndl et ahlliboial f^i. M14 found that the sink 
accretion follows “non-linear stochastic processes”, and con¬ 
cluded that the upper IMF power-law is not a result of BHL 
accretion. 

One of the problems in analyzing the results of numeri¬ 
cal simulations in terms of equation © is that a is not con¬ 
stant, because the medium surrounding accreting objects is 
highly structured and time-variable. In addition, it is not 
trivial to relate the relative velocity in equation o, which 
is presumed to be that of the “ambient medium” far from 
the (isolated) accreting mass, with the motions observed in 
simulations with many gravitating regions, including both 
sinks and gas (see, e.g. iBonnell et aDl2001a^ . Nevertheless, 
if equation © is interpreted as referring to an average over 
volume and time, it may still be relevant for the forma¬ 
tion of massive objects, as gravitational acceleration must 
be present, especially on large scales. 

As numerical simulations of the stellar IMF become 
increasingly sophisticated, with ever-more complex physics 
(e.g., complex equations of state, radiative transfer, stellar 
feedback), it becomes increasingly difficult to isolate the ef¬ 
fects of individual assumptions. These considerations moti¬ 
vated us to take a simplified approach. In this paper, we 
present numerical experiments intended to address the ef¬ 
fects of gravity on the development of upper-mass power-law 
IMFs. We assume isothermality, which allows us to minimize 
the importance of (and uncertainties in) thermal physics, 
and easily produce sufficient numbers of gravitationally- 
bound systems (sink particles) to provide reasonable statis¬ 
tics for the upper-mass IMF slope F. 

Our simulations starting with supersonic but decaying 
turbulence produce an upper-mass sink mass function with 
F = —1 ± O.f; the run starting with density fluctuations 
produced a power-law upper mass tail of comparable slope 
which emerges at late r times. Thus, w e find the power-law 
behavior predicted bv IZinneckerj (Il982l l. even though a sim¬ 
plistic application of equation © with a = constant, does 
not reproduce the evolving nature of the simulations. We 


argue however that these departures from the simple model 
can be reconciled with the simulations by taking into ac¬ 
count the local depletion of matter in the environments of 
the sinks and the mass dependence of the time of sink for¬ 
mation. While our calculations are of a highly simplified 
situation compared with that of real star-forming regions, 
they nevertheless demonstrate the potential of gravitational 
focusing to explain the upper mass power law of the IMF. 


2 METHODS 

2.1 Numerical simulations and initial conditions 


We performed numerical simulations of isothermal gas to 
represent the interior of a small (-^ parsec size) molecular 
cloud. We used the N-bod y, smoothed pa rticle hydrodynam¬ 
ics (SPH) code Gadge t 2 ( Spring^ 20051) wi th the inclusion 
of sink particles as in ijappsen et al.l ll2005ll . Sink particles 
thus represent stars in the sense that they do not have gas 
properties any more, and that the mass that fell into them 
will not emerge back into the medium. Since we do not in¬ 
clude radiative feedback or stellar heating, their interaction 
with the gas was via gravitational forces only. 

Our simulations were performed using 6 million parti¬ 
cles, with a total mass of 1000 Mq, in a cubic box of side 
1 pc in all simulations but one, which was performed in a box 
3 pc on a side. Sink creation was allowed above densities of 
1.7 X 10^cm“®, with an outer radius of Rout = 1.8 x 10“® pc 
within which no further sinks may be made, and an inner 
sink radius of Rin = 1.8 x 10“^ pc or 37 AU. All runs were 
isothermal with an ad opted T = 10 K. According to the 
iBate fc Burkerd d 19971 ) criterion, our mass resolution limit 
is ~ 0.013M©; in practice nearly all of the sink initial masses 
were > O.OOM©. Gas is accreted into sinks if the particles 
are within Rout and are gravitationally-bound to the sink, 
checking previously that another sink does not preferentially 
accrete the particle, or if the gas particle passes within Rtu- 
As our sink particles exert gravity, we impose a gravitational 
smoothing length of 0.003 pc. 

Table 1 shows the main parameters of our runs. Note 
that because these simulation s are isothermal, they can be 
rescaled Isee lHsu et ahlfioiol l. We imposed initial velocity 
fluctuations with a pure rotational veloci ty power spectrum 
with random phases and amplitudes fe.g.. lStone et al]ll998l 'l. 
and with a peak at wavenumber fcfor = Ftt/Lo, where Lq is 
the linear size of the box. No forcing at later times is im¬ 
posed. All calculations but run 33 were evolved for -^1.1 
free-fall times defined by the initial density in the computa¬ 
tional box (see Table 1). 

In the case of run 33, the velocity field was set to zero, 
and initial density fluctuations 5p were imposed, following 
a Kolmogorov Pk oc spectrum with positive ampli¬ 

tudes only. We then stretched these fluctuations such that 
the density field was set to 


_ ^Q(log'5p-(logi5p>)l 


( 2 ) 


Snapshots of the properties of the simulations were 
taken every 10"^ yr, except in the cases of runs 22HR and 32 
for which snapshots were produced every 1000 and 5000 yr, 
respectively. 

We note that all of the runs in Table 1 were actually 
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Table 1. Parameters of the simulations 


Name 

Mach # 

Mass (Mq) 

Size (pc) 

no (cm ®) 

Mj (Mo) 

Tff (Myr) 


Atdump (yr) 

22 

8 

1000 

1 

17,000 

0.48 

0.24 

1.17 

IxlO'^ 

22HR 

8 

1000 

1 

17,000 

0.48 

0.24 

1.08 

1x10® 

23 

16 

1000 

1 

17,000 

0.48 

0.24 

1.00 

1x10"^ 

24 

4 

1000 

1 

17,000 

0.48 

0.24 

1.17 

1x10'^ 

32 

8 

1000 

3 

630 

2.5 

1.26 

1.19 

5x10® 

33 

0 

1000 

1 

17,000 

0.48 

0.24 

0.75 

1x10"^ 


performed twice, with the same parameters but with dif¬ 
ferent initial random seeds. In the following discussion the 
figures generally show the results of one of the two calcula¬ 
tions, except for the mass functions where we add the results 
of both sub-runs to improve the statistics on the slope of the 
mass function. 


2.2 Measuring local properties around the sinks 

Following eq. ©, we will need to calculate the density po 
and velocity dispersion Uturb of the gas around each sink. 
To do so, it is necessary to define some arbitrary radius. In 
principle, it is desirable to calculate these properties outside 
the Bondy-Hoyle radius, given by 

Rbh = 2GM/vtot 

= 8.6 X 10~^ pc (M/Mq) (utot/km sec~^) (3) 

(with Utot = (wturb + Cs) [see eq. ([T], Wturb's the turbulent, or 
non-thermal velocity dispersion of the gas around the star, 
and Ca is the sound speed). As can be seen, while the veloc¬ 
ity dispersion defines the BH radius, it in turn depends on 
the size of the volume at which the measurement is taken. 
Furthermore, since statistically speaking, the velocity dis¬ 
persion increases with size, the larger the volume, the larger 
the velocity dispersion and the smaller the Bondi-Hoyle ra¬ 
dius. To measure Uturb, we have, thus, taken a compromise: 
a volume large enough to ensure to take measurements at 
least at 3 times the BH radius, but small enough to not get 
a BH radius that is arbitrarily small. Thus, we have calcu¬ 
lated the velocity dispersion and the mean density of the 
sph particles that are within I = 6 times the outer radius 
for accretion. Router, and verified that Rbh ^ 3x Router 
in all cases. 


3 RESULTS 

3.1 Initial velocity perturbations 

Figure [T] shows a projection of the evolution of the density 
integrated through the computational box, along with sink 
particles as they form, at four times for run 22. The initial 
supersonic velocity fluctuations shock and rapidly dissipate 
energy, resulting in a (typical) complex distribution of dense 
clumps and filaments. The gas then globally collapses under 
its own self-gravity, the density in the clnmps and filaments 
increase, ultimately forming sink particles which then con¬ 
tinue to accrete from the environment. 

We calculated accretion rates into sinks M as the dif¬ 
ference in sink mass between two snapshots, as a function 



Figure 1. Column density of run 22 at times 0.4, 1.1, 1.8 and 
2.5 xlO® yr. 


of mass. In Figure [2] we show the time sequence of accretion 
rates of sinks as a function of mass (M(sink)). While the 
dependence of the accretion rates on sink mass is crudely 
comparable to the M oc prediction of the simple model, 
there are departures. At the low-mass end of the mass func¬ 
tion (M(sink) ~ O.IMq), there is a large spread in the ac¬ 
cretion rate at a given mass, while at the high mass end, the 
spread in accretion rate vs. mass is much less, but the slope 
is flatter than M~^, with an approximate ® dependence 
at the end of the simulation for the high-mass sinks. These 
departures from the prediction of equation o are due to 
the reduction of matter in the environments of the sinks as 
mass is accreted, and the time dependence of sink creation 
as a function of mass (see discussion in §4). The low-mass 
sinks are particularly affe cted by “competitive accretion” 
(e.g., iBonnell et al.ll^01al f9l. where mass is gravitationally 
attracted to other sinks, particularly the higher-mass ones. 
This competition results in more low-mass sinks exhibiting 
very low accretion rates at later times. 

Figure [3] shows the evolution of the sink mass function 
as a function of time for run 22 (a summation of two in¬ 
dividual runs with different random seed initializations for 
the phases and amplitudes of the velocity fluctuations). In- 


log(N [cm"*]) iog(N [cm'*]) 
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Figure 2. Accretion rates into the sinks M vs. the mass of the 
sinks M(sink) for run 22 from £ = 2.1 to 2.8 x 10^ yr. The dotted 
line indicates a slope dM/dt oc M(sink)^. 


Figure 3. The IMFs for run 22, summed over the two random 
initializations, from t = 2.0 to 2.7x10^ yr (compare with Figure 
[1). The dotted line indicates a power law slope F = — 1. 


terestingly, a quasi-power law distribution of masses is estab¬ 
lished quite early. Most of the sinks are formed at low masses 
~ 0.05 — O.IMq. The slope of the power-law is roughly 
r ~ — 1 over most of the period of sink formation, but it is 
not until the end of the calculation that the statistics become 
good enough to calculate it accurately. Measuring from a low 
mass of O. 2 M 0 to the upper limit yields T = —1.12 ± 0.07 
or from O.SM©, T = —1.15 ± 0.10. 

In Figure |4] we show the final IMFs for the four runs 
23, 24, and 32. are F = -1.24 ± 0.30, -1.04 ± 0.07, and 
—0.96 ± 0.10 for runs 23, 24, and 32, respectively, measur¬ 
ing from the 0.5Mq bin, with essentially the same results 
measuring from 0.2Mq. The first three runs demonstrate 
that changing the Mach number of the initial velocity fluc¬ 
tuations makes essentially no difference to the final results. 
Run 32 consists of the same mass within a volume 27 times 
larger and thus an initial density that factor smaller; the 
resulting IMF is identical within errors to the others. This 
demonstrates the scale-free nature of the problem, with the 
peak in the mass function controlled by the basic resolu¬ 
tion of these isothermal simulations. Note that we achieve 
a F ~ — 1 power-law mass function despite the departures 
of the accretion rates from a strict M oc M{sink)^ relation 

(ED- 

In Figure we use run 32 to explore the onset of the 
sink mass function in more detail. As the initial free-fall time 
of this simulation is a factor of 5.2 longer than in the other 



Moss [Mjyn] 

Figure 4. The final IMFs for individual runs 23, 24, and 32. 

runs, and snapshots were dumped every 5000 yr, this simu¬ 
lation effectively provides higher time resolution (recall that 
the isothermal nature of the simulations permits straight¬ 
forward rescaling). With higher (effective) time resolution, 
the early development of a quasi-power law distribution of 
sink masses is observed clearly. 


3.2 Initial density fluctuations 

Figure shows the surface density for run 33 at t =0.4, 
1.1, 1.8 and 2.5 xl0®yr. The adopted initial conditions are 
such that the density fluctuations are mostly confined to 
a few large-scale concentrations at early times. The ensu¬ 
ing gravitational collapse results in the formation of much 
larger-scale filaments, with much more linear collections of 
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Figure 5. IMFs for run 32 at t = 1.22, 1.3, 1.4 and 1.5 xlO®. 



0.0 0.2 0.4 0,6 0.8 1.0 

X [pc] 


Figure 6. Column density at t = 2.5 X 10®yr for run 33. Black 
dots denote the current positions of the sinks. 

sink particles than in the velocity fluctnation cases (Figure 

m- 

In Figure [3 we show M vs. M{sink) at various times. 
The accretion rates show considerably less correlation with 
sink mass than for the velocity fluctuation case (Figure [2]). 
This difference is probably due to the differing spatial distri¬ 
butions of the sinks, especially at early times. In the density 
fluctuation case, the dense regions initially are far more spa¬ 
tially concentrated into narrow filaments, resulting in higher 
concentration of sinks and thns more competition for acc- 



0.01 0.10 1.00 10.00 0.01 0.10 1.00 10.00 


Moss 

Figure 7. M vs. M(sink) for run 33 from t = 1.1 to 1.8 xlO® 

yr- 

retable mass, causing more scatter in accretion rates at a 
given mass (compare Figure [6] to Figure [T]). 

Figure [8] shows the development of the sink mass func¬ 
tion for run 33. The evolution differs from that of the ve¬ 
locity fluctuation cases in that a quasi-log-normal distribu¬ 
tion emerges initially, consistent with the weaker correlation 
of accretion rates with mass. However, even in this case, a 
roughly power-law tail develops at later times, with an evo¬ 
lution arguably toward F ~ — 1, although the statistics are 
not good enough to make an accurate slope measurement. 

4 DISCUSSION 

4.1 Is BHL accretion relevant? 

The standard BHL accretion equation m does not apply di¬ 
rectly to our simulations because it envisages a well-defined 
external medium from which to accrete which is unper¬ 
turbed by any gravitational field other than that of the 
central mass. This is manifestly not the case because the 
self-gravity of the gas is important and even dominant early 
on, and many other gravitating objects are present. Nev¬ 
ertheless, the development of the power-law mass function 
with r = — 1 as predicted suggests that some basic feature of 
BHL accretion dominates the production of the mass func¬ 
tion even in the more complex situation. 

In the simple BHL model, using equation 0 and as- 
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Figure 8. IMFs for density-fluctuation run 33 from t = 1.1 to 
1.8 xlO® yr. 

suming a = constant, the growth of the accreting mass as a 
function of time is 


where Mo is the mass at t = 0. From this. IZinneckerl (Il982l l 
showed that the resulting mass function Q = dN/dlog M for 
an initial population of masses Mo,i with a mass distribution 
Co = C(^ — 0) is of the form 

* (n^) ^ ® 

More generally, recognizing the environment of the sinks 
is both temporally and spatially variable, one may write the 
accretion rate as 

M = a{t, r)M" . (6) 

Even if a is not constant, if it does not depend on M either 
explicitly or implicitly, t he analysis of th e Boltzmann conti¬ 
nuity equation made bv IZinneckeil (Il982ll is applicable, and 
the asymptotic mass function ( — dN/d\ogM has the form 

0) 

While the mass function slope F = — l.OiO.l that we observe 
in our velocity-fluctuation simulations indicates n = 2, the 
snapshots themselves naively suggest accretion rates for the 
most massive sinks with a shallower mass dependence n ~ 
1.5 (e.g., Figure[2)) which would imply F ~ —0.5. 


It is important to recognize that the mass accretion 
rates of the sinks, integrated over time, must map into the 
mass function. To resolve this apparent contradiction, we 
notice, first of all, that sinks are typically formed in groups: 
in Fig.[9]we show a snapshot of run 22ID2 (upper panel) and 
32 (lower panel) at relatively early times, i.e., when only 64 
and 50 sinks have been formed, respectively {t = 220, 000 
yr, equivalent to 0.91 free-fall times (tff) for run 22ID2, and 
1.325 Myr, equivalent to 1.05 tg for run 32). Thus, at first 
glance it could be arguable that such groups share similar 
environmental properties (density and velocity dispersion). 

We also note that, since a = a{t,r), local variations in 
a may blur any quadratic dependence of M on M at a given 
time. Thus, by dividing the mass accretion rate of each sink 
by its local a, it might be possible to recover any power-law 
dependence (eq. [6]) of the accretion mass rate with mass, at 
least by groups with similar environmental properties (i.e., 
similar a): 

log M = log Qi (t, r) + n log M, (8) 

In Fig. [To] we now plot the raticQ M/a' against the mass M 
for all the sinks at the same times shown in Fig. (21 In both 
cases we notice that there are two main groups, each one 
sharing the approximate tendency M oc M^, but with an 
offset of ~3 orders of magnitude between them. This result 
shows that different groups of sinks might have different en¬ 
vironmental properties, and thus, any possible relationship 
of the kind M oc might be blured when the groups are 
analyzed all together. 

To worsen the case, we furthermore notice that even 
when one can account for such kind of segregation by ai, the 
situation might still not be straightforward. For instance, in 
the upper panel of Fig. llOl lwhich corresponds to run 22ID2) 
we denote by empty squares all the sinks at f = 1.03tff 
but those that belong to the group that it is zoomed in in 
the upper panel of Fig. The 16 sinks that belong to that 
group are denoted also with cyan asterisks if the value of 
a' is smaller than and by red crosses otherwise. We 

then notice in the upper panel of Fig. 1101 that all but one 
sinks that belong to the mentioned group have a' < 10^®, 
and that only one of them (red symbol) has M/a' > 

This sink is also denoted by a cross in the zoom of the upper 
panel of Fig.|9](the remaining ones are denoted by asterisks), 
and even though it is the neighborhood of the group, in the 
upper panel of Fig. llOl it falls far away from the M/a' oc M^ 
tendency delineated by the rest of the group. The reason 
might be either that the local density of this sink is smaller 
than that of the rest of the group (e.g., it is located in the 
periphery of the core), or that somehow this sink has a larger 
velocity (e.g., the sink its been ejected), compared to the 
rest of the group, increasing the velocity dispersion of the 
sph particles of its neighborhood. In fact, a careful analysis 
of the zoom in the upper panel of Fig. |9]shows that it might 
be the first case: the sink is located in the periphery of the 
core, according to the column density map. 

In a similar way, we notice that each one of the four 
groups of sinks shown in the lower panel of Fig. (9] also ex¬ 
hibit the approximate M oc tendency, although with 
some scatter, (see lower panel of Fig. dOj), and even with 

^ In what follows, since o. oc p/vtot^ we will define a.' = p/ntot^ 
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Figure 9. Maps of runs 22ID2 (a, upper panel) and 32 {(h)^ 
lower panel) at t = 220,000 yr and at 1.325 Myr, respectively. 
In both cases sinks are formed by groups, where apparently, the 
local density and velocity dispersion of the gas is nearly the same 
for all the sinks in the group. Closer examination shows that this 
is not necessarily the case (see text). 


some outliers (e.g., blue cross in group 1 showing M/a' ^ 
5 — 8x10“^°). These results suggest that even within a sin¬ 
gle group, there migh be different values of Utot and p, and 
then, strong scatter in the plane (M, M/a') migh be seen. 

To further stress this point, in Fig. we plot again 
M/a' vs. M for the same two runs, but (a) at later times, 
when substantially more sinks have been formed (363 for run 
22ID2, upper panel, and 354 for run 32, lower panel), and 
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Figure 10. M/{pv^) vs. M for a) run 22ID2, upper panel) and 
b run 32 at the same times shown in Fig. [9] The dashed line has 
a slope of 2. Note that when one accounts for the local variations 
of a', even though there is some scatter, individual groups ex¬ 
hibit nearly the quadratic dependence of M with M. The groups 
labeled in the lower panel are the same than those shown in the 
lower panel of Fig. |9] 


then the density and velocity helds are even more complex, 
and (b) in both cases, we color and label every point accord¬ 
ing to the local value of a' (the actual values are colored 
and labeled in log scale in the inner box). Although there 
is some scatter at every a', it is clear that in each case we 
recover the quadratic dependency per groups, strongly sug¬ 
gesting that, indeed, the accretion is of Bondi-Hoyle type, 
with non-constant a, as would be predicted by eq. (|H1) for 
n = 2 and different a. In summary, even for a single group, 
substantial scatter in the plane (M, M/a') can be due to 
local fluctuations in velocity or density. These fluctuations 
blur the quadratic dependence of M/a' with M. 

Since individual groups show the quadratic dependence 
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Figure 11. M/a' vs. M for runs a 22ID2 at f = l.OStff and b 
32 t = 1.17tfr. The dashed line has a slope of 2, and the actual 
values of o' are labeled and colored accordingly in the rectangle, 
in units of [cm km sec“^]~®. Note that when one accounts for 
the local variations of o', individual groups exhibit the quadratic 
dependence of M with M. 


when we take the ratio Mja, why does the plot M vs. M 
for all the groups exhibits a shallower slope? Should not 
they all collapse to a single M oc relation, maybe with 
substantial scatter, when plotting M vs. Ml As shown in 
Fig. 1121 where we plot the final mass of each sink vs. the 
time at which every sink is formed, the most massive sinks 
systematically form earlier than the lowest-mass sinks, and 
thus accrete for longer. In other words, the integration of M 


10.000 

^ '.000 
2 “ 

U) 

g 0.100 
2 

0.010 

0.001 

1.1 1.2 ',3 1.4 1.5 

time of creation [Myr] 

Figure 12. Time at which sinks initially form for Run 32. The 
sinks with larger final masses systematically accrete for longer 
than the lowest-mass sinks (see text). Similar plots are found for 
all the runs. 




time [Myr] 

Figure 13. Mass growth vs, time for sinks in run 32, separated 
into final mass bins for clarity. Upper left, sinks with eventual 
masses > IMq; upper right, sinks with final masses 1.0 < M ^ 
O. 5 M 0 ; lower left, 0.1 < M ^ O.SMq; lower right, M ^ 0.1 sinks. 


over time has an implicit dependence on mass. This allows 
the mass function to evolve to F = — 1. Thus, as time pro¬ 
ceeds, the local environment becomes depleted as it goes into 
sinks, with the result that earlier-formed sinks accrete some¬ 
what more slowly than when they were created. This envi¬ 
ronmental depletion is translated in the shallower slopes of 
individual sink mass evolution vs. time (see Figure [13} than 
predicted by equation (|3} with a' =cs10. To furth er support 
our co nclusion, we note that the simulations of IHsu et al.l 
(|201(]|) . where all sinks were formed at the same time and 
continuing sink creation was not allowed, do show M oc 
for the massive sinks. 

The self-gravitating nature of the medium, unlike that 


^ Note, however, that by imposing a decreasing a{t) as a function 
of time, one can find individual mass accretion histories similar 
to those shown in Fig. 1131 
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Figure 14. Zoom-in of the time evolution of a region in run 22 
which produces one of the earliest, and ultimately most massive, 
sinks. The sink (when formed) is indicated by the large symbol at 
(0,0). The red particles are those that will be accreted during the 
next time step. The “hole” around the sink corresponds to the 
region R < R(inner), inside of which it is assumed all gas par¬ 
ticles accrete. Gravitationally-driven infall via streams is evident 
(see text). 


envisaged in the simple BHL calculation, is responsible for 
the typically rapid mass growth at early times, as in the “col¬ 
lapse” mode identified by M14. It is not surprising that this 
is qualitatively different than that predicted by equation ©, 
as this early evolution represents the initial formation of the 
gravitating body, and, as such, is not addressed by the BHL 
formula, which assumes the presence of an already existing 
mass. One must also recognize that the initial sink mass 
is only the central condensation of a larger self-gravitating 
clump. In essence, the early evolution is that of protostel- 
lar core collapse; the initial protostellar core has a much 
lower mass than its surrounding envelope, but the latter is 
bound by its own self-gravity, not that of the protostellar 
core. The lack of applicability of gravitational focusing at 
early stages, however, does not mean that gravitational fo¬ 
cusing or quasi-BHL accretion cannot occur later to build 
the mass well beyond initial collapse values. 

One may ask how this differs from models in which tur¬ 
bulence produces an initial clump ma ss that sets the final 
mass , even for massive objects (e . g., ^doan fc NordlunJ 
I 2 OO 2 I ; [Hennebelle fc Chabrier|[200a . |2009I L We illustrate the 
difference in Figure 1141 where we show the evolution of a 
region around a sink that forms early and becomes massive. 
In the upper left panel, the red particles are those that will 
become part of the initial sink in the next time step (the 
sink has not yet formed). As described in the previous para¬ 
graph, the initial sink mass is only a fraction of the larger 
self-gravitating clump. As time proceeds, the added mat¬ 


ter arises from gravitationally-driven streams which do not 
represent initial turbulent fragments. 

We note in passing that even though the “velocity fluc¬ 
tuation” models were initiated with a supersonic turbulent 
velocity field, the motions rapidly decay, leaving density 
fluctuations as the dominant features which then produce 
gravitationally-induced supersonic motions. Thus, it may be 
that initial clump masses are also influenced by gravitational 
focusing. 


4.2 Comparison with other investigations 

As noted in the Introduction, in contrast to our hndings, 
M14 concluded that “the upper power law tail of the IMF 
is unrelated to Bondi-Hoyle accretion”, and, instead, argued 
that the mass function is the result of a distribution of (ini¬ 
tial) seed masses, growth times, and stochastic, fluctuating 
accretion rates with M oc We suggest that the pri¬ 

mary reason for the differing conclusions is that their mass 
function is considerably different from ours. The mass func¬ 
tion of M14 (their Figure 11) shows a much broader peak 
than ours; inspection of their Figure 7 shows that most of 
the objects in the broad peak of the mass function corre¬ 
spond to “collapse-dominated” sinks, with the mass accreted 
in only one or two episodes. These objects arguably corre¬ 
spond to “initial” fragmentation masses, and so BHL accre¬ 
tion is undoubtedly irrelevant for these objects (see previous 
section). However, we also note that the M14 mass func¬ 
tion at masses above IMq shows an approximate power-law 
slope of r ~ —1.35. M14 label these “strong” accretors or 
“accretion-dominated” systems, and we suggest that BHL- 
type accretion is actually occuring in this mass range. 

The broad peak in the M14 mass function relative to 
ours may reflect, in part, the pie c ewise barotropic equa¬ 
tion of state used bv iBonnell et al. l (|2008l). which can affect 
fragmentation (see discussion in IJaPDsen et"^ (l2005l B. Our 
adoption of an isothermal equation of state means that the 
peaks of our low-mass mass functions are essentially set by 
our resolution limits rather than by any thermal physics; this 
likely expands the dynamic range in mass over which BHL 
accretion can operate, resulting in a clearer development of 
the Zinnecker power law in our calculations. 

As argued in the previous section, the apparent weaker 
dependence of accretion rate on mass than oc M(sinfe)^ does 
not rule out BHL accretion because of the depletion of en¬ 
vironmental gas along with competition between sinks. In 
fact, when accounting for local variations of a, the quadratic 
dependence emerges. While M14 obtain an even slower de¬ 
pendence of mass accretion on sink mass than we do, this is 
partly the result of anchoring the ht at masses below IMq 
( their Figure 6), where accretion beyond the initial fragmen¬ 
tation mass is small. 

Finally, M14 also called attention to rapid, large fluc¬ 
tuations in mass accretion rates as a possible problem with 
BHL accretion; such variations have been seen in other simu - 
lations as well (IPeters et al.ll20f^ . l201ll : |Padoan et al.l[2014|j . 
We also see strong variations in M as well, but this variabil¬ 
ity does not prevent the development of the F = — 1 power 
law. These accretion rate variations reflect clumping in the 
near-sink environments (see Figure 03 . It appears that these 
clumps are regions which are attempting to fragment and 
gravitationally collapse, but which fail to do so before being 
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accreted. In our runs, additional sink formation near exist¬ 
ing sinks is prohibited to avoid numerical problems and short 
time steps as new sinks begin to orbit each other closely. As 
a test, we reran one of the simulations for a short period of 
time, removing this constraint, and indeed found additional 
sinks would form. Thus, it is not clear whether the “inges¬ 
tion” of this material really represents accretion of gas that 
will be accreted into the central object, or really represents 
fragmentation into an additional sink that will simply end 
up orbiting the first sink. (Of course there is a general is¬ 
sue in this type of simulation, in that resolution limits mean 
that we are only reporting “system” or “enclosed masses 
within some radius”.) In our case, the assumption of isother- 
mality probably enhances the tendency toward gravitational 
clumping. Thus, even though highly variable accretion may 
be of great interest for under standing the time-v ariability 
of protostellar accretion (e.g. lAudard et al.ll2014h . further 
careful consideration of missing effects such as stellar feed¬ 
back, magnetic fields, etc. is needed before applying these 
results to real systems. 


5 CONCLUSIONS 

Our numerical experiments with simplihed physics indi¬ 
cate that gravitational focusing is an effective and plausible 
means for producing an upper mass power law T ~ — 1. The 
situation is more co mplicated t han in the simplistic analysis 
of BHL accretion bv IZinneckeil lll982fl , but even the temporal 
and spatial variability, and the presence of nearby gravitat¬ 
ing masses, does not prevent the development of the upper 
mass power law. Although the details of accretion vs. time 
and sink mass are complicated, we argue that the relative 
accretion rates, intergrated over time and space, exhibit an 
dependence as in BHL accretion, which produces the 
power-law mass function we find. 

The present set of simulations are obviously quite ide¬ 
alized. Our purpose was to elucidate the basic physics of 
gravitationally-focused accretion in as easily-visualized and 
interpretable a situation as possible. While effects of mag¬ 
netic fields, non-isothermal behavior, and stellar feedback 
are clearly important to any realistic and comprehensive un¬ 
derstanding of the stellar mass function, the results shown 
here strongly suggest that BHL accretion, understood in its 
most general sense as an “attractor” with an dependence 
on relative accretion rates, remains relevant to an under¬ 
standing o f the upper mass regio n of the IMF, as prev iously 
argued by iBonnell et al.l (l2001al [9i ; IClark et al.l (l2007li . 
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